---
title: "How Civil Resistance Improves Inclusive Democratization: Replication of Main Paper Analyses"
author: "Subindra Bogati, Titik Firawati, Jonathan Pinckney, Ches Thurber"
date: "April 2024"
output:
  html_document: default
bibliography: grateful-refs.bib
---

This document produces all tables and figures appearing in the main text of "How Civil Resistance Improves Inclusive Democratization," as prepared for final publication in *Perspectives on Politics* in April 2024. It requires the set of packages listed in the code as well as the data file "IncDemData_April2024.csv". Raw data sources and the code used to create that data file can be obtained by the authors.

The tables and figures can be replicated by knitting the .Rmd file to .html output. Alternatively, each "chunk" can be run independently in RStudio. The output as published was produced using R version 4.3.2. Additional software packages are listed below.

Please contact Ches Thurber (ches.thurber@gmail.com) with any questions regarding data and replication.

## Software Used
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)

library(dplyr) #v1.1.4
library(tidyr) #v1.3.1
library(ggplot2) #v3.5.1
library(ggeffects) #v1.3.0
library(knitr) #v1.43
library(kableExtra) #v1.4.0
library(modelsummary) #v1.4.2
library(vdemdata) #v14.0
library(sensemakr) #v0.1.4
library(grateful) #v0.2.4
```

```{r, Cite packages, message = F, warning = F}
pkgs <- cite_packages(output = "table", out.dir = ".")
```

|Package      |Version |Citation                                       |
|:------------|:-------|:----------------------------------------------|
|base         |4.3.2   |@base                                          |
|ggeffects    |1.3.0   |@ggeffects                                     |
|kableExtra   |1.4.0   |@kableExtra                                    |
|knitr        |1.43    |@knitr2014; @knitr2015; @knitr2023             |
|modelsummary |1.4.2   |@modelsummary                                  |
|plm          |2.6.3   |@plm2008; @plm2017; @plm2018                   |
|rmarkdown    |2.24    |@rmarkdown2018; @rmarkdown2020; @rmarkdown2023 |
|sensemakr    |0.1.4   |@sensemakr                                     |
|tidyverse    |2.0.0   |@tidyverse                                     |
|vdemdata     |14.0    |@vdemdata                                      |

## Figure 1: Trends in Annual Averages of VDem Inclusion Measures, 1945-2014

```{r, inclusion}
vdem <- vdemdata::vdem
vdem_inclusion <- vdem %>% 
  select(COWcode, year, v2x_egal, v2xpe_exlgender, v2xpe_exlecon, v2xpe_exlsocgr) %>%
  filter(year>1944 & year<2015)

Inc_Trends <- vdem_inclusion %>%
  group_by(year) %>%
  summarize(
    Equality = mean(v2x_egal, na.rm=T),
    Gender = 1-mean(v2xpe_exlgender, na.rm = T),
    Class = 1-mean(v2xpe_exlecon, na.rm = T),
    Ethnicity =  1-mean(v2xpe_exlsocgr, na.rm = T)
  ) 

trends_long <- Inc_Trends %>%
  gather("id", "value", 2:5) %>%
  mutate(Type = factor(id, levels = c("Equality", "Gender","Class", "Ethnicity")))


myColors <- c("#FF0000", "#1B9E77", "#7570B3", "#E6AB02" )
names(myColors) <- c("Equality", "Class","Gender", "Ethnicity")

myLines <- c("solid", "twodash", "longdash", "dotdash")

ggplot(trends_long, aes(year, value, color = Type)) +
         geom_line(aes(linetype = Type)) +
  labs(y = "Inclusion Index", x = "Year", color = "Inclusion Measure", linetype = "Inclusion Measure") +
  scale_x_continuous(limits = c(1940, 2020), labels = c(1940, 1960, 1980, 2000, 2020)) +
  theme_bw() +
  scale_color_manual(values = myColors) 

ggsave("Figure1.png", width = 6, height = 3, units = "in")

```

## Table 1: Average Changes in Inclusion Measures One Year and Five Years after Transition

```{r, desc stats}

IncDemData <- read.csv("IncDemData_April2024.csv")

Table1 <- IncDemData %>%
  group_by(crt) %>%
  summarize(
    Total = n(),
    Equality_DY1 = mean(equality_d1, na.rm=T),
    Equality_DY5 = mean(equality_d5, na.rm=T),
    Gender_DY1 = mean(gender_d1, na.rm=T),
    Gender_DY5 = mean(gender_d5, na.rm=T),
    Class_DY1 = mean(class_d1, na.rm=T),
    Class_DY5 = mean(class_d5, na.rm=T),
    Ethnicity_DY1 = mean(ethnicity_d1, na.rm=T),
    Ethnicity_DY5 = mean(ethnicity_d5, na.rm=T)
  )

kbl(Table1, booktabs = T, digits = 3) 
#%>%
  #kable_styling(latex_options = "scale_down")
```



## Table 2: OLS Regressions of Effect of Transition Type on Changes in Inclusion, Before vs. One Year After Transition

```{r,  1-year regressions}

Model1 <- lm(equality_d1 ~ crt  + duration + polyarchy_l1 + gdppc_l1  + civsoc_l1 + violence + prev_CRT + CW + equality_l1 + equality_ld5, data = IncDemData)

Model2 <- lm(gender_d1~crt  + duration + polyarchy_l1 + gdppc_l1 + civsoc_l1 + violence + prev_CRT + CW + gender_l1 + gender_ld5, data = IncDemData)

Model3 <- lm(class_d1~crt  + duration + polyarchy_l1 + gdppc_l1  + civsoc_l1 + violence + prev_CRT + CW + class_l1 + class_ld5 , data = IncDemData)

Model4 <- lm(ethnicity_d1~crt  + duration + polyarchy_l1 + gdppc_l1  + civsoc_l1 + violence + prev_CRT + CW + ethnicity_l1 + ethnicity_ld5, data = IncDemData)


Table2 <- list(
                  "Model 1: Equality" = Model1,
                  "Model 2: Gender" = Model2,
                  "Model 3: Class" = Model3,
                  "Model 4: Ethnicity" = Model4)

modelsummary(Table2, stars = T, coef_map = c("crt" = "CR Transition", "equality_l1" = "Prior Inclusion", "gender_l1" = "Prior Inclusion", "class_l1" = "Prior Inclusion", "ethnicity_l1" = "Prior Inclusion", "equality_ld5" = "Prior Trend", "gender_ld5" = "Prior Trend", "class_ld5" = "Prior Trend", "ethnicity_ld5" = "Prior Trend", "duration" = "Duration", "polyarchy_l1" = "VDem Polyarchy", "civsoc_l1"="Civil Society", "gdppc_l1" = "GDP pc (log)", "violence" = "Violence", "prev_CRT" = "Prior CRT", "CW" = "Cold War"), gof_map = c("nobs", "r.squared"), vcov = "stata", notes = "HC1 robust standard errors used across models. Prior Inclusion and Prior Trend covariates are specific to the inclusion measure DV and vary between models. CR Transition is a dichotomous variable indicating a Civil Resistance Transition (CRT) as compared to a non-Civil Resistance Transition (non-CRT). ")

```

## Figure 2: Coefficient estimates (marginal effects) from full models
```{r, Fig 2 Coefficients}
modelplot(list("Ethnicity" = Model4, "Class" = Model3, "Gender" = Model2 , "Equality" = Model1),  coef_map = c('crt'='CR Transition'), vcov="stata") +
  aes(shape = model, color = model) +
    labs(x = 'Coefficients', 
         y = 'Covariates',
         color = "Inclusion Measure",
         shape = "Inclusion Measure") +
  scale_x_continuous(breaks=c(0, 0.02, 0.04, 0.06, 0.08, 0.10)) +
  geom_vline(xintercept=0) +
  theme_bw() +
  scale_shape_manual(values = c(15, 21, 17, 19), guide = guide_legend(reverse = TRUE)) +
   scale_color_manual(values = myColors, guide = guide_legend(reverse = TRUE))
  
ggsave("Figure2.png", width = 6, height = 3, units = "in")
```

## Figure 3: Expected changes in inclusion measures for CRTs versus non-CRTs
```{r EV plot}
Model1_EV<- ggeffect(Model1, terms=c("crt"), vcov.fun = "vcovHC", vcov.type = "HC1")
Model1_EV$Dimension <- c("Equality", "Equality")
Model1_EV$x <- c("Non-CRT", "CRT")

Model2_EV<- ggeffect(Model2, terms=c("crt"), vcov.fun = "vcovHC", vcov.type = "HC1")
Model2_EV$Dimension <- c("Gender", "Gender")
Model2_EV$x <- c("Non-CRT", "CRT")

Model3_EV<- ggeffect(Model3, terms=c("crt"), vcov.fun = "vcovHC", vcov.type = "HC1")
Model3_EV$Dimension <- c("Class", "Class")
Model3_EV$x <- c("Non-CRT", "CRT")

Model4_EV<- ggeffect(Model4, terms=c("crt"), vcov.fun = "vcovHC", vcov.type = "HC1")
Model4_EV$Dimension <- c("Ethnicity", "Ethnicity")
Model4_EV$x <- c("Non-CRT", "CRT")

Total_EV <- rbind(Model2_EV, Model3_EV, Model4_EV, Model1_EV)
Total_EV$Dimension <- factor(Total_EV$Dimension, levels = c("Equality", "Gender","Class","Ethnicity"))

ggplot(Total_EV, aes(x, predicted, color=Dimension)) +
  aes(shape = Dimension) +
    geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
   labs(x = "Transition Type",
         y = "Expected Change in Inclusion",
        color = "Inclusion Measure",
        shape = "Inclusion Measure") +
  facet_wrap(~Dimension) +
    theme_bw() + 
  scale_shape_manual(values = c(19, 17, 21, 15)) +
  scale_color_manual(values = myColors)

ggsave("Figure3.png", width = 6, height = 3, units = "in")
```

## Table 3: OLS Regressions of Effect of Transition Type on Changes in Women’s Political Representation and Mass Mobilization

```{r,  Mechanism Tests}

Model5 <- lm(gen_polpart_d1 ~ crt  + duration + polyarchy_l1 + gdppc_l1  + civsoc_l1 + violence + prev_CRT + CW + gen_polpart_l1 + gen_polpart_ld5, data = IncDemData)


Model6 <- lm(mobilization_d1~crt  + duration + polyarchy_l1 + gdppc_l1  + civsoc_l1 + violence + prev_CRT + CW + mobilization_l5, data = IncDemData)



Table3 <- list(
                  "Model 5: Women's Political Participation" = Model5,
                  "Model 6: Post-Transtion Mass Mobilization" = Model6
                  )
                 

modelsummary(Table3, stars = T, coef_map = c("crt" = "CR Transition", "gen_polpart_l1" = "Prior Leg. Rep.","gen_polpart_ld5" = "Prior Leg. Rep. Trend" , "mobilization_l5" = "Prior (5yr) Mobilization", "duration" = "Duration", "polyarchy_l1" = "VDem Polyarchy", "civsoc_l1"="Civil Society", "gdppc_l1" = "GDP pc (log)", "violence" = "Violence", "prev_CRT" = "Prior CRT", "CW" = "Cold War"), gof_map = c("nobs", "r.squared"), vcov = "stata", notes = "HC1 robust standard errors used across models.")


```



## Table 4: Average Changes in Inclusion Measures by Group Participation


```{r, Table 4, warning = F , message= F}

gen_div_sum <- IncDemData %>%
  filter(crt==1) %>%
  group_by(div_gender) %>%
  summarize(Gender_N = n(),
            Gender_DY1 = mean(gender_d1, na.rm=T),
    Gender_DY5 = mean(gender_d5, na.rm=T))%>%
  mutate(Type = ifelse(div_gender == 0, "CRT No Participation", "CRT w/ Participation"))

class_div_sum <- IncDemData %>%
  filter(crt==1) %>%
  group_by(div_class) %>%
  summarize(Class_N = n(),
            Class_DY1 = mean(class_d1, na.rm=T),
    Class_DY5 = mean(class_d5, na.rm=T)) %>%
  mutate(Type = ifelse(div_class == 0, "CRT No Participation", "CRT w/ Participation"))

eth_part_sum <- IncDemData %>%
  filter(crt==1) %>%
  group_by(excl_part) %>%
  summarize(Ethnicity_N = n(),
            Ethnicity_DY1 = mean(ethnicity_d1, na.rm=T),
    Ethnicity_DY5 = mean(ethnicity_d5, na.rm=T)
            ) %>%
  mutate(Type = ifelse(excl_part == 0, "CRT No Participation", "CRT w/ Participation"))

Table4 <- left_join(gen_div_sum, class_div_sum, by = "Type") %>%
  left_join(eth_part_sum, by = "Type") %>%
  select(Type, Gender_N:Gender_DY5, Class_N:Class_DY5, Ethnicity_N:Ethnicity_DY5)

kbl(Table4, booktabs = T, digits = 3) 

```



## Table 5: Regressions of Campaign Participation and Transition Type

```{r, warning = F, message = F}

IncDemData$gen_type <- factor(IncDemData$gen_type)
Model9 <- lm(gender_d1~gen_type + duration + polyarchy_l1 + gdppc_l1 + civsoc_l1 + violence + prev_CRT + CW + gender_l1 + gender_ld5, data = IncDemData)

IncDemData$class_type <- factor(IncDemData$class_type)
Model10 <- lm(class_d1~class_type + duration + polyarchy_l1 + gdppc_l1 +  civsoc_l1 + violence + prev_CRT + CW + class_l1 + class_ld5, data = IncDemData)

IncDemData$eth_type <- factor(IncDemData$eth_type)
Model11 <- lm(ethnicity_d1~eth_type + duration + polyarchy_l1 + gdppc_l1  + civsoc_l1 + violence + prev_CRT + CW +  ethnicity_l1 + ethnicity_ld5, data = IncDemData)

Table5Models <- list(
                  "Model 9: Gender" = Model9,
                  "Model 10: Class" = Model10,
                  "Model 11: Ethnicity" = Model11)

modelsummary(Table5Models, stars = T, coef_map = c("eth_type1" = "CRT No Participation", "eth_type2" = "CRT w/ Participation",  "class_type1" = "CRT No Participation", "class_type2" = "CRT w/ Participation" , "gen_type1" = "CRT No Participation" , "gen_type2" = "CRT w/ Participation", "gender_l1" = "Prior Inclusion", "class_l1" = "Prior Inclusion", "ethnicity_l1" = "Prior Inclusion", "gender_ld5" = "Prior Trend", "class_ld5" = "Prior Trend", "ethnicity_ld5" = "Prior Trend", "polyarchy_l1" = "VDem Polyarchy", "civsoc_l1"="Civil Society", "gdppc_l1" = "GDP pc (log)", "violence" = "Violence", "prev_CRT" = "Prior CRT", "CW" = "Cold War"), gof_map = c("nobs", "r.squared"), vcov = "stata", notes = "HC1 robust standard errors used across models. Prior Inclusion and Prior Trend covariates are specific to the inclusion measure DV and vary between models." )
```

## Figure 4: Pairwise marginal effects of contrasts

```{r, contrasts}
IncDemData$gen_type_reorder <- factor(IncDemData$gen_type, levels = c(2,1,0))
GenPart_Reorder<- lm(gender_d1~gen_type_reorder + duration + polyarchy_l1 + gdppc_l1 + civsoc_l1 + violence + prev_CRT + CW + gender_l1 + gender_ld5, data = IncDemData)
PC_Gen <- hypothesis_test(GenPart_Reorder, terms = "gen_type_reorder")
PC_Gen<- PC_Gen %>%
  mutate(Dimension = c("Gender", "Gender", "Gender"),
         Pair = c("CRT: No Part -> CRT: Yes Part.", "Non-CRT -> CRT: Yes Part", "Non-CRT -> CRT: No Part")) %>%
  select(-gen_type_reorder)

IncDemData$class_type_reorder <- factor(IncDemData$class_type, levels = c(2,1,0))
ClassPart_Reorder <- lm(class_d1~class_type_reorder + duration + polyarchy_l1 + gdppc_l1 +  civsoc_l1 + violence + prev_CRT + CW + class_l1 + class_ld5, data = IncDemData)
PC_Class <- hypothesis_test(ClassPart_Reorder, "class_type_reorder")
PC_Class <- PC_Class %>%
  mutate(Dimension = c("Class", "Class", "Class"),
         Pair = c("CRT: No Part -> CRT: Yes Part.","Non-CRT -> CRT: Yes Part", "Non-CRT -> CRT: No Part"))%>%
  select(-class_type_reorder)
           
IncDemData$eth_type_reorder <- factor(IncDemData$eth_type, levels = c(2,1,0))
EthPart_Reorder <- lm(ethnicity_d1~eth_type_reorder + duration + polyarchy_l1 + gdppc_l1  + civsoc_l1 + violence + prev_CRT + CW +  ethnicity_l1 + ethnicity_ld5, data = IncDemData)
PC_Eth <- hypothesis_test(EthPart_Reorder, "eth_type_reorder")
PC_Eth <- PC_Eth %>%
  mutate(Dimension =  c("Ethnicity", "Ethnicity", "Ethnicity"),
         Pair = c("CRT: No Part -> CRT: Yes Part.","Non-CRT -> CRT: Yes Part", "Non-CRT -> CRT: No Part")) %>%
  select(-eth_type_reorder)

PC_Full <- rbind(PC_Gen, PC_Class, PC_Eth)
PC_Full$Dimension <- factor(PC_Full$Dimension, levels = c( "Ethnicity","Class","Gender"))

ggplot(PC_Full, aes(Contrast, Pair, color=Dimension, shape = Dimension)) +
    geom_pointrange(aes(xmin = conf.low, xmax = conf.high), position = position_dodge(width = 0.4)) +
  geom_vline(xintercept = 0) +
    labs("Pairwise Marginal Effect",
         y = "Contrast",
         color = "Inclusion Measure",
         shape = "Inclusion Measure") +
  theme_bw() +
  scale_shape_manual(values = c(15, 21, 17, 19), guide = guide_legend(reverse = TRUE)) +
   scale_color_manual(values = myColors, guide = guide_legend(reverse = TRUE))

ggsave("Figure4.png", width = 6, height = 3, units = "in")

```

## Figure 5: Expected changes in inclusion measures for CRTs with excluded group participation versus without versus non-CRTs

```{r Figure 5}


GenPartFull1_EV<- ggeffect(Model9, terms=c("gen_type"), vcov.fun = "vcovHC", vcov.type = "HC1")
GenPartFull1_EV$dimension <- c("Gender", "Gender", "Gender")
GenPartFull1_EV$x <- c("No CRT", "CRT: No Part.", "CRT: Yes Part.")

ClassPartFull1_EV<- ggeffect(Model10, terms=c("class_type"), vcov.fun = "vcovHC", vcov.type = "HC1")
ClassPartFull1_EV$dimension <- c("Class", "Class", "Class")
ClassPartFull1_EV$x <- c("No CRT", "CRT: No Part.", "CRT: Yes Part.")

EthPartFull1_EV<- ggeffect(Model11, terms=c("eth_type"), vcov.fun = "vcovHC", vcov.type = "HC1")
EthPartFull1_EV$dimension <- c("Ethnicity", "Ethnicity", "Ethnicity")
EthPartFull1_EV$x <- c("No CRT", "CRT: No Part.", "CRT: Yes Part.")

PartFull1_EV <- rbind(GenPartFull1_EV, ClassPartFull1_EV, EthPartFull1_EV)
PartFull1_EV$dimension <- factor(PartFull1_EV$dimension, levels = c("Equality", "Gender","Class","Ethnicity"))

ggplot(PartFull1_EV, aes(x, predicted, color=dimension, shape = dimension)) +
    geom_pointrange(aes(ymin = conf.low, ymax = conf.high)) +
   labs(x = "Transition Type",
         y = "Expected Change in Inclusion") +
  facet_wrap(~dimension) +
    theme_bw() +
  scale_shape_manual(values = c(15, 21, 17, 19), guide = "none") +
   scale_color_manual(values = myColors, guide = "none")
ggsave("Figure5.png", width = 10, height = 5, units = "in")
```




## Figure 6: Changes to T-value for Civil Resistance Transition coefficient estimate from a hypothetical omitted confounder

```{r, Sensitivity}
sensitivity <- sensemakr(model = Model1,
                         treatment = "crt",
                         benchmark_covariates = "equality_l1",
                         kd=1:3)

summary(sensitivity)
```


```{r, Figure 6}

plot(sensitivity, sensitivity.of = "t-value")

png(file = "Figure6.png",
width=6, height=6, units = "in", res = 1200)
plot(sensitivity, sensitivity.of = "t-value", xlim = c(0,0.25), ylim = c(0,0.25))
```

## Complete Software Package Citations